function out = loopMatRegr(data,settings)

returns = data.returns;
spreads = data.spreads;
mats = data.mats;
rec = data.rec;
q = data.q; 
h = data.h;

NWlags = 2*h;
blockSize = NWlags;

doBoot = settings.doBoot;
Nsim = settings.Nsim;
recLb = settings.recLb;

idxnonan = ~isnan(rec); 

if ~isempty(data.q);
      
   idxnonan = idxnonan & ~isnan(q);   
   q = q(idxnonan);  
   
end

rec = rec(idxnonan);
spreads = spreads(idxnonan,:);
returns = returns(idxnonan,:);


T = size(returns,1);

numk = length(mats);

for j = 1:numk %loop over maturities

     %Newey-West std errors
     [bv, sebv, r2,~,fit,out.residuals(j,:)] = robustOLS(returns(:,j),[ones(T,1) spreads(:,j) q],NWlags,1);

     out.mat_bv(j,:) = bv';
     out.mat_sebv(j,:) = sebv';
     out.mat_r2(j,:) = r2;
     out.mat_fit(j,:) = fit;

     e_j = returns(:,j)-[ones(T,1) spreads(:,j) q]*bv;
     e0_j = returns(:,j)-mean(returns(:,j));
     
     out.mat_r2rec(j,:) = 1-(e_j(rec==1)'*e_j(rec==1)) / (e0_j(rec==1)'*e0_j(rec==1));
     out.mat_r2exp(j,:) = 1-(e_j(rec==0)'*e_j(rec==0)) / (e0_j(rec==0)'*e0_j(rec==0));
     
     [bvSwitch,sebvSwitch,r2Switch,~,fitSwitch] = robustOLS(returns(:,j),[1-rec spreads(:,j).*(1-rec) rec spreads(:,j).*rec q],NWlags,1);

     out.mat_bvSwitch(j,:) = bvSwitch';
     out.mat_sebvSwitch(j,:) = sebvSwitch';
     out.mat_r2Switch(j,:) = r2Switch;
     out.mat_fitSwitch(j,:) = fitSwitch;

     [bvDiff,sebvDiff,r2Diff] = robustOLS(returns(:,j),[ones(T,1) spreads(:,j) rec spreads(:,j).*rec q],NWlags,1); 

     out.mat_bvDiff(j,:) = bvDiff';
     out.mat_sebvDiff(j,:) = sebvDiff';
     out.mat_r2Diff(j,:) = r2Diff;

     eDiff_j = returns(:,j)-[ones(T,1) spreads(:,j) rec spreads(:,j).*rec q]*bvDiff;
     
     bv0Diff = robustOLS(returns(:,j),[ones(T,1) rec],NWlags,1); 
     e0Diff_j = returns(:,j)-[ones(T,1) rec]*bv0Diff;

     out.mat_r2recDiff(j,:) = 1-(eDiff_j(rec==1)'*eDiff_j(rec==1)) / (e0Diff_j(rec==1)'*e0Diff_j(rec==1));
     out.mat_r2expDiff(j,:) = 1-(eDiff_j(rec==0)'*eDiff_j(rec==0)) / (e0Diff_j(rec==0)'*e0Diff_j(rec==0));         
     
     out.mat_tv(j,:) = bv'./sebv';
     out.mat_tvDiff(j,:) = bvDiff'./sebvDiff';    
     out.mat_tvSwitch(j,:) = bvSwitch'./sebvSwitch'; 
         
     if doBoot == 1
     %Stationary bootstrap
     [bvBoot, sebvBoot] = statBootMinRec( returns(:,j), [ones(T,1) spreads(:,j) q], ...
                                                                             rec, recLb, Nsim, blockSize );

     out.mat_sebvBoot(j,:) = sebvBoot';
     out.blockSize = blockSize;

     [bvSwitchBoot, sebvSwitchBoot] = statBootMinRec( returns(:,j), [1-rec spreads(:,j).*(1-rec) rec spreads(:,j).*rec q], ...
                                                                                  rec, recLb, Nsim, blockSize );

     out.mat_sebvSwitchBoot(j,:) = sebvSwitchBoot';         

     [bvDiffBoot, sebvDiffBoot] = statBootMinRec( returns(:,j), [ones(T,1) spreads(:,j) rec spreads(:,j).*rec q], ...
                                                                            rec, recLb, Nsim, blockSize );

     out.mat_sebvDiffBoot(j,:) = sebvDiffBoot';  
     
     out.mat_tvBoot(j,:) = bv'./sebvBoot';
     out.mat_tvDiffBoot(j,:) = bvDiff'./sebvDiffBoot';   
     out.mat_tvSwitchBoot(j,:) = bvSwitch'./sebvSwitchBoot';   
          
     
     end
     
     
end



end